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ABSTRACT 



Context. Runaway massive stars are O- and B-type stars with high spatial velocities with respect to the interstellar 
medium. These stars can produce bowshocks in the surrounding gas. Bowshocks develop as arc-shaped structures, with 
bows pointing to the same direction as the stellar velocity, while the star moves supersonically through the interstellar 
gas. The piled-up shocked matter emits thermal radiation and a population of locally accelerated relativistic particles 
is expected to produce non-thermal emission over a wide range of energies. 
Aims. We aim to model the non-thermal radiation produced in these sources. 

Metliods. Under some assumptions, we computed the non-thermal emission produced by the relativistic particles and 
the thermal radiation caused by free-free interactions, for 04I and 091 stars. We applied our model to ( Oph (HD 
149757), an intensively studied massive star seen from the northern hemisphere. This star has spectral type 09.5V and 
is a well-known runaway. 

Results. Spectral energy distributions of massive runaways are predicted for the whole electromagnetic spectrum. 
Conclusions. We conclude that the non-thermal radiation might be detectable at various energy bands for relatively 
nearby runaway stars, especially at high-energy gamma rays. Inverse Compton scattering with photons from the heated 
dust gives the most important contribution to the high-energy spectrum. This emission approaches Ferrm sensitivities 
in the case of Oph. 

Key words, stars: early-type - gamma-rays: theory - radiation mechanisms: non-thermal - stars individual: ^ Oph. 
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1. Introduction 

Runaway stars have high peculiar velocities, v^, ^ 30 km 
s-i (e.g. Gies & Bolton 1986; Tetzlaff, Neuhauser & Hohle 
2011). These stars move supersonically through the inter- 
stellar medium (ISM) and are ejected from their birth as- 
sociations by gravitational kicks (e.g. Fujii & Zwartz 2011). 
Runaway stars with strong winds - spectral types OB and 
Wolf-Rayet (WR) ~ can produce stellar bowshocks when 
interacting with the surrounding gas (e.g. Van Buren et al. 
1995). Bowshocks develop as arc-shaped structures, with 
bows pointing ahead of the stars (in the same direction as 
the stellar velocity). 

There are several cataloged objects of this type (e.g. Van 
Buren et al. 1995; Noriega-Crespo, van Buren, & Dgani 
1997; Kobulnicky, Gilbert & Kiminki 2009; Gvaramadze, 
Kniazev, & Kroupa 2011; Peri et al. 2011). Not ah run- 
away stars produce a bowshock. Van Buren et al. (1995) 
concluded that only 30% of the cataloged runaway stars 
present bowshock-like structures. More recently. Peri et al. 
(2011) found out that only ^ 10% of early-type runaway 
stars develop an observable bowshock in the infrared (IR) . 
Some physical situations exist in which a bowshock does 
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not form. If the star is moving in a low-density, hot plasma 
then it might be moving subsonically and a shock will not 
occur. Also, a bowshock might not form if the star has a 
weak wind or if it has a high space velocity (see Comeron 
& Kaper 1998, Huthoff & Kaper 2002). 

Bowshocks occur around many classes of astrophysical 
sources: pulsars, cataclysmic variables, colliding wind bina- 
ries, cometary H II regions, and even in groups and clusters 
of galaxies. A lot of work has been done on bowshock mod- 
eling (e.g. Van Buren & McCray 1988; Van Buren, Mac 
Low, Wood & ChurchweU 1990; Bandiera 1993; Van Buren 
1993; Brighenti & Dercole 1995; Wilkin 1996; Comeron 
1997; Chen & Huang 1997; Comeron & Kaper 1998; Wilkin 
2000; Wareing, Zijlstra & O'Brien 2007). Related physi- 
cal problems such as solving the astrophysical blunt-body 
problem are treated in e.g. Canto & Raga (1998) and more 
recently in Schulreich & Breitschwerdt (2011); the physics 
of mixing layers is discussed in e.g. Baranov, Krasnobaev & 
Ruderman (1976), Canto & Raga (1991) and Raga, Cabrit 
& Canto (1995). 

The supersonic stellar wind sweeps the ISM material 
and piles it up in the bowshock. The stellar and shock- 
excited radiation heats this swept-up material. The dust, 
in turn, re-radiates the energy as mid-to-far IR flux (Van 
Buren & McCray 1988). 

Relativistic particles can be accelerated at strong shocks 
producing non-thermal emission (Drury 1983). Benaglia et 
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Fig. 1. Scheme of the system of shocks (not to scale). The 
acceleration region is indicated by dots. 



al. (2010) have reported non-thermal radio emission from 
the bowshock of the runaway star BD +43° 3654. This emis- 
sion is thought to be synchrotron radiation generated by the 
interaction of relativistic electrons with the magnetic field 
of the source. 

In this work we present a model for the radiative non- 
thermal emission that takes place in the bowshocks of run- 
away stars. We compute the emission produced by the rel- 
ativistic particles accelerated at the shock, and the thermal 
emission caused by free- free mechanism (Bremsstrahlung). 
Finally, we apply this model to the well-known stellar bow- 
shock from the 09.5V star ( Oph. The bowshock shape is 
also computed, following the analytical method developed 
by Wilkin (2000). 

Our paper is organized as follows. In the next section we 
present the radiative model developed to compute the non- 
thermal emission and the absorption that can take place 
in these types of sources. We also present a brief discus- 
sion on particle acceleration in shocks. In Sec. |3] we ap- 
ply the model to a particular source, the bowshock from C 
Oph. We present the best-fit shell shape and the computed 
spectral energy distribution. Finally, in Sec. |4] we discuss 
the detectability of the non-thermal radiation and offer our 
conclusions. 



2. Radiative model 

Bowshocks from runaway stars produce thermal emission 
by radiative heating of the swept-up dust by the stellar ra- 
diation field. The infrared signal of the heated dusty bow- 
shock is strongest in the far infrared (e.g. Van Buren & 
McCray 1988; Kobulnicky, Gilbert & Kiminki 2009). The 
dust temperature can be understood in terms of dust mod- 
els. For a typical bowshock of a runaway star, a fraction ~ 
10~^ of the star bolometric luminosity is emitted in the in- 
frared. In this paper we do not calculate this radiation, and 
leave it for a future work (del Valle et al., in preparation). 

The shocked ISM can also produce thermal emission 
through free- free interactions (Bremsstrahlung). This emis- 
sion peaks at energies ~ 1 eV; we calculate it in Sec l2.5.2l 

The interactions of locally accelerated relativistic par- 
ticles with the matter, radiation and magnetic fields in the 
system produce non-thermal radiation. We calculate this 
emission below. 



2.1. Shocks and particle acceleration 

Shocks transfer kinetic energy to non-thermal particles 
through particle acceleration. The acceleration mechanism 
is diffusive shock acceleration, the so-called Fermi I mech- 



anism (e.g. Bell 1978). In this mechanism particles are ac- 
celerated by successive bouncing across the shock, gaining 
energy in each cross. The particle defiection is mediated 
by magnetic field irregularities. These irregularities are ex- 
pected from turbulent and magnetic instabilities. For the 
mechanism to operate is necessary that the particles in the 
shocked medium (downstream) can effectively diffuse and 
reach the shock. 

The particle energy gain in each cycle (upstream- 
downstream- upstream) is AE/E o: (ws/c), where Vs is the 
shock velocity; and after k cycles the particle energy is 

E — Ei(^l + ^) , where Ei is the initial energy. The time 
that a particle requires to accelerate up to an energy E is 
given by 



eBc 



(1) 



Here, B is the magnetic field in the acceleration region, and 
T] is the acceleration efficiency (Drury 1983): 



IQ D \ c) 



(2) 



D is the diffusion coefficient, and rg — E/{eB) is the parti- 
cle gyroradius. In the Bohm limit Db = rgc/3. Fast shocks 
are more efficient accelerators than slow shocks. The result- 
ing spectrum of injected particles by this mechanism is a 
power-law, i.e. Q{E) oc (e.g. Protheroe 1999), with 

a -2- 2.2. 

The collision of the supersonic stellar wind with the ISM 
around a runaway star results in a system of two shocks 
(see Fig. [1]). In the steady state, mass and momentum are 
conserved and a ffow is established between the two shocks 
that carries away the mass and momentum deposited by the 
colliding fluids. In the reference frame of the star, the ISM 
can be considered as an incoming wind of parallel stream- 
lines and the stellar wind as a radial outflow (e.g. Wilkin 
1996). One of the shocks is a forward shock that propagates 
in the same direction as the stellar wind with velocity V[s 
~ Vi,; and the other is a reverse shock that propagates in 
the opposite direction with Vrs ~ Vw. The stellar wind can 
be considered as a continuous power source, therefore both 
shocks reach a steady state. If radiative losses are ineffi- 
cient within the shock discontinuity, the shock is adiabatic; 
otherwise, the shock is radiative. 

The stellar wind is divergent because its ram pressure 
decreases with distance, while the ram pressure of the ISM 
is constant. The point where the ram pressure of the wind 
and the ISM balance, i.e. PwV^ = PaV^, is called the stand- 
off point, where pw — M^/AttR^Vv,. This point of the shell 
defines the standoff radius Rq: 



4^PaK' ■ 



(3) 



The physical conditions in bowshocks from runaway 
stars produce an adiabatic and fast reverse shock at the 
shocked stellar wind; and the forward shock is radiative 
and slow (e.g. Van Buren 1993). The radiative shock com- 
presses the material, and as the temperature decreases, the 
density grows. 

Kis et al. (2004) demonstrated that upstream particles 
in the Earth bowshock undergo diffusive transport into the 
upstream region, a direct evidence of Fermi I acceleration 
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(e.g. Burgess 2007). Benaglia et al. (2010) reported the de- 
tection of non-thermal radio emission from the bowshock of 
the runaway star BD-H43°3654. The non-thermal radiation 
is expected from synchrotron emission generated by rela- 
tivistic electrons accelerated at either the forward shock or 
the reverse shock. 

We considered that an initial supra-thermal population 
of relativistic particles, electrons and protons, are acceler- 
ated in the adiabatic shock (reverse shock). This shock is 
faster than the forward shock, so in the reverse shock the 
relativistic particles are accelerated more efficiently by the 
Fermi mechanism (see Eq. All equations given above 
are valid for planar shocks. We considered an acceleration 
region where the shock is nearly flat (see FiglT]). The width 
of the shocked stellar wind. A, can be estimated as A ^ 
M^^Rq, where M is the Mach number of the shocked wind. 
The acceleration region is assumed as a small region near 
the bowshock apex, of scale length ~ A; the acceleration 
volume is given in Table [TJ 

The available power in the system is the kinetic power 
from the stellar wind: 



(4) 



To estimate the magnetic field in the flow, we consider 
that the magnetic energy density is in subequipartition with 
respect to the kinetic energy Lt, by a 0.1 factor, i.e.: 

(5) 



Stt 



O.ILt 



where A is the area of a sphere of radius Rq. This condition 
ensures that the flow is matter-dominated, i.e. compress- 
ible, allowing shocks to develop. 

The kinetic power available in the acceleration region is 
L = /Lt, where / is the ratio of the volume of a sphere of 
radius Rq and the volume of the acceleration region. A small 
fraction of this kinetic power goes into relativistic particles, 
-^rei = QreiL. We adopted a standard fraction ^rei = 0.1 (e.g. 
Protheroe 1999). We took into account both hadronic and 
leptonic content in the relativistic power, L^^i = Lp + aL^. 
The ratio of relativistic protons to electrons, a, is unknown. 
We considered two cases a = 1 (equal energy density in 
both species) and a = 100 (as observed in Galactic cosmic 
rays, Ginzburg & Syrovatskii 1964). 

Here we discuss two different types of massive star: an 
041 and an 091 star. The values for the parameters involved 
are given in Table [1] 



2.2. Non-thermal radiative losses 

The electrons lose energy mainly by inverse Compton 
(IC) scattering, synchrotron radiation, and relativistic 
Bremsstrahlung. 

The synchrotron cooling time rate is 



tz 



4 ffTcU B f TTl, 
3 TTlpC^ 



E 



m J mc'^ 



(6) 



where ctt is the Thomson cross section and Ub is the mag- 
netic energy density. 

The IC losses can be calculated from (Blumenthal & 
Gould 1970) 



(7) 



where e and ei are the incident photon and scattered photon 
energy, respectively, and 



dA^ 1 27rrQ mc^ nph(e)d6 



dMei Ec 7 
with 

f{q)^2q\nq+{l 



2.)(l-.) + i^(l 
2 a + bq 



9), 



(8) 



(9) 



where b = Aej/mc^ and q = ei/[b{Ea — ei)]. This expres- 
sion takes into account losses in the Klein-Nishina regimen. 
The target photon fields are the radiation fields of the ac- 
celeration region. The fields considered here are the stellar 
radiation field at a distance Rq from the star - assumed 
as a black body at temperature T^,- , and the IR radiation 
from the heated dust, also considered as a black body at 
TiR. Dust grains are heated by starlight, and cool by radi- 
ating in the infrared. To estimate Tir we used a simplified 
dust model based on Draine & Lee (1984): 



rTR = 27 a-y^xye ^-1/3 k 



*38 ^'■Opc 



(10) 



(e.g. Van Buren & McCray 1988). Here a^m ~ 0.2 fim is the 
dust grain radius. More detailed and complex dust emission 
models are beyond the scope of this work. For more details 
see e.g. Draine & Li (2007) and Draine (2011). 

The relativistic Bremsstrahlung losses are calculated 
considering a complete ionized plasma, using (Berezinskii 
et al. 1990) 



In 



2E, 

771c 



1 



(11) 



where 7i is the density of target ions in the acceleration 
region (the shocked stellar wind in this case). The density 
of the shocked wind according to the Rankine-Hugoniot 
equations for adiabatic shocks (e.g. Landau & Lifshitz 1959) 
is 47i„, where ri^ is the wind density. 

Protons lose energy through proton-proton inelastic col- 
lisions with the shocked wind material. These interactions 
produce neutral and charged pions; the former decay and 
produced gamma rays (e.g. Vila & Aharonian 2010). The 
latter decay into secondary electrons/positrons and neu- 
trinos. The losses produced by proton-proton interactions 
occur on time scales of 



jicdppKpp , 



(12) 



where n is the density of target protons, ~ 4nw, and Kpp is 
the ineslasticity (^ 0.5). The cross section can be approxi- 
mated by (Kelner, Aharonian, & Bugayov 2006) 



app = (34.3 + 1.88L + 0.25L^ 



1 - 



th 



Er, 



(13) 



where L = ln(i;p/lTeV); and Eth = 280 MeV. 

Photomeson production is irrelevant at the energies con- 
sidered in this paper. 

In Fig. 12.21 we show the cooling rates for electrons and 
protons in the acceleration region for an 041 and an 091 
star. In the 041 star system the IC scattering of IR photons 
dominates the radiative losses. For the 091 star, the IC 
scattering of IR photons and synchrotron radiation prevail 
among the radiative losses. 
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Parameter 


04 


09 


Ro Standoff radius [pc] 


8.3 


0.2 


Vt Spatial velocity [km s~^] 


100 


30 


Wind velocity [km s"^] 


2.2x10^ 


0.8x10^ 


Wind mass loss rate [Mq yr""*^] 


10~* 


lO"*^ 


/Za Ambient density [cm""^] 


1 


100 


R'' Macrnetir field [HI 


^ Q n V 1 




Tj Acceleration efficiency 


~ z.U X iU 


~ z. / X iU 


L Available power [erg s 


~ O.Z X lU 


~ 4.0 X lU 


n rTa n Ton-1"ri- Ipt^I^oti pnPT'Cf^r r^\'\c\ 


I 


100 


Qrei Jet content of relativistic particles 


10% 


10% 


a Injection index 


2 


2 


A Thickness of shocked wind [Ro] 


~ 0.3 


~ 0.3 


Volacc Acceleration region volume [cm^"^] 


~ TxlO^** 


~ lO'^i 


L*'^ Star luminosity [Lq] 


~ 7x10^ 


~ 5x10" 


Tib. Dust temperature [K] 


~ 24 


~ 54 



Table 1. Parameters for the different types of stars. 
'^Values derived by Kobulnicky, Gilbert & Kiminki (2009). 

^This value corresponds to the magnetic field in the acceleration region, obtained from Eq. ([5]). 
'^Values from Martins, Schaerer & Hillier (2005). 



The relativistic particles can also suffer from non- 
radiative losses due to escape from the acceleration region. 
Particles can be convected away by the stellar wind on a 
time tconv ~ ^/Vw This non-radiative loss dominates the 
proton cooling rate in both types of runaway stars (041 and 
091). 

The minimum kinetic energy for each particle is con- 
sidered to be on the order of the rest mass energy. The 
maximum energy for the electrons and protons is obtained 
equating the smallest cooling rate to the acceleration rate, 
given by Eq. ([TJ. In the 041 system the electrons reach 
energies ^10 TeV and the protons can reach energies ^ 
10^ TeV. For 091 electrons reach energies ~ TeV, and pro- 
tons energies ^ 10 TeV. In both cases the Hillas criterion 
is satisfied: 



i^max < 300 (rg/cm) (B/G) eV, 



(14) 



The normalization constant Qq for each type of particles 



is obtained from L^.p as 



(17) 



Here V is the volume of the acceleration region. 

The exact solution to the equation is also a power-law: 



N{E) = 



with 



dE 
"dF 

xexp 



dE'Q{E') 



loss"^ E 

t{E,E') 



r dE" 


dE" 


IE 


dt 



(18) 



(19) 



loss 



where rg is the maximum particle gyroradius in the avail- 
able space, i.e. ~ A. 

2.3. Particle distributions 

To calculate the steady state particle distributions N{E) 
for electrons and protons, we solved the transport equation 
in steady state (Ginzburg & Syrovatskii 1964): 



d_ 

'dE 



dE 
'At 



N(E) 



N{E) 



to. 



Q{E), 



(15) 



where ^csc is the convection time. We assumed that all phys- 
ical properties in the acceleration region are homogeneous. 

Q{E) is the injection function, a power-law in the energy 
of the particles, as expected from diffusive shock accelera- 
tion: 



Q{E)^QoE-". 



(16) 



2.4. Non-tliermal radiative processes 

As mentioned before, the interactions of the relativistic par- 
ticles with the various fields in the system produce non- 
thermal emission. 

We considered synchrotron emission, inverse Compton 
scattering with the stellar photon field and the IR emis- 
sion, and relativistic Bremsstrahlung from electrons. For 
protons, inelastic collisions between the relativistic protons 
with the shocked wind material were also calculated. 

The synchrotron emission is given by 

n-Emax 



KssA{E.y)E^V 



1.85 X \^ 



hmc^ 

1/3 



En 



dEN{E)^ 
Ec 



(E^ 



where 



3 ehB ( E 



Air mc 



(20) 



(21) 
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Fig. 2. Acceleration and cooling time scales for electrons and protons for an 041 star (up), and for an 091 star (down). 
Left panels are for electrons and right panels are for protons. The non-radiative losess produced by convection are also 
shown. 



'^SSA is the synchrotron self-absorption (SSA) factor: 

1 _ g— tssaC-E^) 



KSSA 



TSSa(^^7) 



(22) 



here tssa(£'^) is the SSA optical depth (see Rybicki 
Lightman 1979). 

The IC emission is calculated using 



dePic{E,,E^,e). 

The spectrum of scattered photons is 

p /p. p. N 3crTc(meC^)^ nph(e) 
Pic{Ec,E^,e) = — — f{q), 

with f{q) given by Eq. ([9]). 

The relativistic Bremsstrahlung luminosity is 

L^{E.,) = E^V / naBiE„E^)—N,{E,)dE,, 
Je-, 47r 

where 

aB{Ee,E^) = -^cl>{E,,E^), 



(23) 



(24) 



(25) 



(26) 



and 

(t>{E,,E^) = [1 + (1 - E^/E,f ~ 2/3(1 - E^/E,)] 



x<\n 



2E,{E, -E^) 1 
mcC^E^ 2 



(27) 



To compute the gamma-ray emission produced by neu- 
tral pion decays for Ep < 0.1 TeV, we used the following 
expression: 



L^{E^) = 2VE^ [ 



q^riE-^) 



^E^ - mlc^ 



dE^, 



(28) 



with i?min — E^ + nriT^c^ / AE-„. The Tr'^-emissivity is given 
by (Aharonian & Atoyan 2000) 

QTviE^) = ^^pp (^rripC^ + Et,/k.^ Jp ^TOpC^ -I- E^/kJ^ (29) 

with - 0.17 (Gaisser 1990). For Ep < 0.1 TeV down to 
the threshold energy, the following replacement is necessary 

6{E K^-Ekin). (30) 

Here h is the total number of tt'^ created per p — p collision, 
and i^kin = Ep — nipC^ is the proton kinetic energy. 

The gamma-ray luminosity in the range 0.1 TeV< Ep < 
10^ TeV can be obtained from (Kelner et al. 2006): 

/>oo 

L^{E^) = nE^^^V / (J,,,,x[Ep)Np{Ep) 
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(31) 



with Fjy-^,Epj a function of and Ep (see Vila & 

Aharonian 2009 and references therein). 

Secondary electron-positron pairs are created through 
the p—p inelastic collisions (see Orellana et al. 2007). These 
particles can also produce non-thermal emission, which we 
took into account. For the relevant formulae see Orellana 
et al.'s paper. 

Figure [3] shows the computed SED for both types of 
stars. We calculated the contribution of the secondary pairs, 
to synchrotron emission and IC luminosity, only for the 091 
star (because we assumed a proton dominated content in 
this case). The non-thermal luminosity of the 041 star is 
higher, mainly because the available power in the system is 
greater. The IC emission produced in the interaction with 
the dust photons dominates both SEDs at high energies. 
At lower energies synchrotron radiation dominates, but is 
relatively more important in the 091 system. 

The IC emission from the dust photons, in general, in- 
creases as the temperature of the dust also increases. This 
temperature depends, as a first aproximation, on and Rq 
(see Eq. ([TU])). As argued in Sec. 12.11 B must be at sube- 
quipartition values. This parameter not only affects the syn- 
chrotron emission, but the acceleration rate (see Eq. ((!])). 
A higher value of B gives a lower acceleration rate and the 
maximum particle energies change (decrease). 

2.5. Emission of the shocked ISM 
2.5.1. Non-thermal emission 

Some relativistic particles can diffuse from the accelera- 
tion region to the shocked ISM before losing most of their 
energy. The diffusion time can be estimated as A^/2Db, 
where Db is the diffusion coefficient in the Bohm limit as 
before. Particles that diffuse to the denser shocked ISM (see 
Fig. [1} can interact there through p — p inelastic collisions 
and relativistic Bremsstrahlung. The shocked ISM is con- 
tained within a very thin layer. 

Because of the radiative losses the temperature de- 
creases and the density increases in the shocked ISM; both 
quantities are described by a temperature profile Tfs{x) and 
a density profile nfs{x). The temperature profile is given by 
(Zhekov & Palla 2007) 



-7 X 10" 



-19 



2P. 



ad 



(3.6)a; + T' 



ad 



1/(3.6) 




Log(£) [eV] 

Fig. 4. 7— ray emission from the shocked ISM produced 
hy p — p interactions from an 091 star. 

function injection (see Aharonian & Atoyan 1996; Bosch- 
Ramon, Aharonian & Paredes 2004): 



Q{E,R,t)^Np6{R)6{t), 



(34) 



where iVp is the initial energy distribution of the particles 
accelerated at the reverse shock. The general solution is 



iVpdiff(^;) 



KE-f' exp-E/E„ 
X exp ' 



(35) 



wher 



'pp 



6 X 10^(ri/lcm ■^) ^ yr. Here we consider an 
average value for the shocked density fifg = J^" nfs{x)dx/xc; 
t is the propagation time after the injection into the shocked 
ISM, t ~ cXc- RdiS is the so-called diffusion radius. For t 



<< Tpp, RdiS = 2\/Dt, with D the diffusion coefficient, we 
considered D = Db, see Sec. 12.11 

The particle energy distributions of the diffused elec- 
trons can be approximated by 



Afcdiff-Q(S,i?,t)icooi, (36) 

where 1/icooi = l/tr/nchr + *Tg + tsLmss- The interac- 
tions of diffused particles with the matter of the shocked 
ISM -through p — p interactions for protons and relativistic 
Bremsstrahlung for electrons- produce non-thermal emis- 
sion. 

For the two stars considered here, only the protons from 
the 091 system can diffuse into the shocked ISM. Fig. S] 
shows this contribution to the non-thermal SED. This emis- 
sion peaks at E ^ 100 MeV, with L ^ 10^^ erg s~^. From 
Fig. [3] it can be seen that this contribution is negligible 

(32) 

compared to the total IC radiation. 



where x is the distance from the shock; Tad ^ 2 x 10^ K 
and Pad = jPa^l, are the corresponding values of temper- 
ature and pressure for the adiabatic case. 
The density profile is 



nisix) 



_ C'o^ad 

kBT{s{x) ' 



(33) 



with Co a normalization constant such that = nts(Ta), 
where Ta is the ambient temperature. We estimated the 
width of the shocked region as Xc, where T{xc) = T^- 

The particle energy distribution of the diffused pro- 
tons is obtained solving Eq. ([T5|) . considering a 3-D delta 



2.5.2. Thermal emission 

The shocked ISM produces thermal emission through free- 
free mechanism (thermal Bremsstrahlung). For complete- 
ness we computed this contribution to the total SED. 

We calculated the thermal emission by integrating the 
emissivity e along the shocked ISM, considering the tem- 
perature and density profiles given by Eq. (15^ and ([55]) . 
The emissivity is given by e.g. (Lang 1999) 

e - 5.4 X 10"3^^^5(i/,r)exp-/izy/fcBr 



(37) 



T' 



ergs ^cm Hz ^ rad ^, 
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Fig. 3. Non-thermal emission from the shocked stellar wind for an 041 star (left), and for an 091 star (right). 




Log(£) [eV] 

Fig. 5. Thermal emission from the shocked ISM produced 
by free-free mechanism for an 041 star and for an 091 star. 

where §{1/, T) is the free-free Gaunt factor given by ~ 
0.54 In [5 X 10^ (T'^/^/i^)]; m and rie are the ion and elec- 
tron densities respectively, we assumed n,; — 12^= nis- The 
emission volume is a spherical wedge of width Xc and ra- 
dius A, see Fig. [TJ In Fig. [5] the result of this contribution 
is shown. 

For the 041 star the free-free emission is negligible (see 
Fig. [3]). For the 091 case, at energies around ~ 1 eV this 
thermal contribution dominates over synchrotron radiation. 

2.6. Absorption 

Gamma rays can be absorbed in the acceleration region 
immediately after they are created, by photon-photon an- 
nihilation. All radiation fields in the acceleration region, 
thermal and non-thermal, provide target photons for pair 
creation. 

The differential opacity for a gamma ray traveling in the 
direction due to photons of an energy e in the direction 
Cph is 



dTj.y = (1 — e^eph)neCr-y^dedr2dZ, 



(38) 



where dfl is the solid angle of the emitting surface and is 
the radiation density. The photon annihilation cross-section 
is (Gould & Schreder 1967) 



r(/3) 



2 



W - 2) 



(3 



(39) 



where /3 = (1 - 1/s)^/^, and s — E^e^l — e^epi^) / (rrieC^)'^ . 
Here, Ej and e are the energies of the gamma ray and the 
target photon, respectively. The threshold energy is given 
by 



2(mec^)^ 

£/^e = r 

(1 - e-yGph) 

The opacity is 



(40) 



(1 



r(/3)nph(e)duded/.(41) 



Here, u = cos??, is the angle between the momenta of 
the colliding photons, and I is the photon path across the 
target radiation field (here 1 — 6). The target photon fields 
are those generated within the acceleration region, the IR 
emission from the heated dust, and the star radiation field 
at a distance d = Rq — S. 

The optical depth is a trajectory integral on which the 
angular dependence has a very significant effect, and conse- 
quently the absorption depends strongly on the line of sight 
(e.g. Romero, del Valle, & Orellana 2010). Photons travel- 
ing from the acceleration region toward an observer placed 
at A, see Fig. [51 are not absorbed by the stellar photon 
field. On the other hand, photons traveling toward an ob- 
server placed at B strongly interact with the stellar photon 
field. This interaction depends on the distance d between 
the star and the photon path, and therefore depends on the 
angle a. 

Photons of lower energies are absorbed by matter 
through photoionization. The target material can be the 
shocked ISM and the material along the line of sight corre- 
sponding to each particular source. The optical depth t^h 
can be approximated as 



yH ^ NHcrjj<i{Ej). 



(42) 



Here Nn is the column density corresponding to the target 
density. The cross section ct^n was taken from Ryter (1996). 
For more details see Reynoso, Medina & Romero (2011) and 
references therein. 

Below we apply the radiative model described so far to 
the specific case of ( Oph. 
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Fig. 6. Diagram of the path followed by a gamma ray trav- 
eling from the acceleration region toward observers placed 
in front and behind the runaway star (not to scale). 



3. Application to ( Oph 

The star C, Oph (HD 149757) is one of the brightest mas- 
sive stars in the northern hemisphere and has been inten- 
sively studied. This star has spectral type 09.5V and is a 
well-known runaway star. It is rapidly rotating with almost 
break-up velocity, with v sin(i) ~ 400 — 500 km s^^ (Walker 
et al. 1979; Repolust et al. 2004). 

C Oph bowshock has been observed by the Infrared 
Astronomical Satellite IRAS (Van Buren & McCray 1988, 
Noriega- Crespo et al. 1997) and more recently by the Wide- 
field Infrared Survey Explorer WIS^ The WISE image 
shows a very clear, regular structure. This source is quite 
nearby, located at ~ 222 pc from the Earth (Megier et al. 
2009). 

3.1. Bowshock shape 

We adopted a mass loss rate for C Oph of ^ lO^*" Mq 
yr~^ (FuUerton et al. 2006) and terminal wind velocity 

= 1550 km s-i (e.g. Hurbig, Oskinova & Scholler 2011). 
For the ISM density we adopted a value of ria ~ 10 cm^'^; 
since ^ Oph is embedded in a cirrus cloud region, the den- 
sity and ambient temperature might be higher than average 
(e.g. Vidal et al. 2011). These parameters yield Rq ~ 0.3 
pc, which agrees well with the upper limit measured by 
Peri et al. (2011). The list of the values of the main param- 
eters adopted in our calculations (we follow Marcolino et 
al. 2009) are in Tabled 

To compute the bowshock shape we used the analytical 
method developed by Wilkin (2000). We assumed that the 
angular dependence of the wind momentum flux is equato- 
rial (C2 = -1 and A = 30°, see Sec. 4.3 of Wilkin 2000). 
Fig. [7] shows the 3-D computed bowshock shape. 

To compare the 2-D image observed by WISE we ro- 
tated the 3-D coordinate system by three angles. We de- 
fined the coordinates Xp = E and Zp = N (see Fig. [5]) to 
describe the WISE image. As usual, the origin is placed at 
the position of the star. From the WISE image of the bow- 
shock of ^ Oph it is evident that the midpoint from the star 
to the bowshock (see Fig. 4 from Peri et al. 2011) is in the 
direction of the star velocity. 

Figure [8] shows the best fit of ( Oph bowshock. 

^ http : I /wise2.ipac.caltech.edu/ docs/release/prelim/ 



?[pcl 




Fig. 7. C, Oph computed bowshock. 
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Fig. 8. Top: image of C Oph from WI SE 
(http : / /wise2.ipac.caltech.edu/docs/release/prelim/). 
Dots pointed by arrows correspond to the projection to 
the 3-D structure. Bottom: ^ Oph projected bowshock. 



Our result shows that a simple model can provide a good 
representation of the observational data, assuming that the 
IR image is a good tracer of the actual hydrodynamic bow- 
shock shape. The differences between the observed and the 
intrinsic shape of the bowshock depend on many factors 
involving radiative transfer, cooling time, dust characteris- 
tics, dust distribution, and so on. 
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Table 2. Parameters for ( Oph 
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Fig. 9. Electron and proton losses, acceleration rates, dif- 
fusion and convection time scales - defined in Sec. 12.51 - for 

COph. 




Fig. 10. Internal opacity. 



3.2. Spectral energy distribution 

In Fig.|n]we show the radiative losses, the acceleration rates 
and the diffusion and convection times (see Sec. 12. 5p for 
electrons and protons. The maximum energy is ~ TeV for 
both species of particles. These values are in accordance 
with the Hillas criterion (see Sec. 12. ip . 

The internal photon-photon optical depth in the bow- 
shock is shown in Fig llOl It is negligible in the energy ranges 
of interest. The photoelectrical absorption is also negligible 
due to the small amount of material that photons cross on 
their way to the observer. The external absorption is also 
negligible given the relative positions of the bowshock, the 
star and the observer. 

Figure [Til shows the computed spectral energy distribu- 
tion (SED) for the emission from the bowshock of ( Oph, 
along with the sensitivity of the gamma-ray detectors CTA 
(Cherenkov Telescope Array-forthcoming-), MAGIC and 
Fermi, the X-ray satellite XMM-Newton (theoretical upper 
limit from Hasinger et al. 2001), and VLA (upper limit from 
the NVSS survey - Condon, Cotton, Greisen et al. 1998 -, 
angular resolution is given by Peri et al. 2011). For com- 



Parameter 



value 





Standoff radius 


0.3 pc 


Mw 


Wind mass loss rate 


10-^ Mq yr-i 


a 


Hadron-to-lepton energy ratio 


1 




Content of relativistic particles 


10% 


a 


Particle injection index 


2 




Wind velocity 


1.5x 10* cm s"^ 


L 


Available power 


5x10^^ erg s"^ 


B 


Magnetic field 


5x10-* G 


V* 


Star velocity 


30 km s"^ 




ISM number density 


10 cm"^ 




Star temperature 


3.2x10'' K 


R* 


Star radius 


9 Rq 




Star luminosity 


10^ Lq 


Tm 


Dust temperature 


~ 66 K 



pleteness the IR IRAS data are also shown (Van Buren & 
MacCray 1988). 

The expected non-thermal luminosity of the source is 
weak. However, since C, Oph is very nearby, the bowshock 
might be detectable at gamma-ray and X-ray wavelengths 
through long exposures, under the assumptions we made. 
We remark that the sensitivity shown in Fig.[TT]for MAGIC 
is for 50 hours of exposure over the source, and for Fermi 
it is for one year of integration. For these types of sources 
an instrument like CTA gives the best chance of detection. 
CTA might become a unique tool to explore the high-energy 
radiation produced by runaway massive stars and the pop- 
ulation of relativistic particles generated in them. 



4. Discussion and conclusions 

Bowshocks of runaway massive stars are natural candidates 
for particle acceleration. The different types of massive stars 
have different energetics (depending on the wind param- 
eters). The available power for particle acceleration also 
depends on the distance between the star and the accel- 
eration region. The different available powers produce dif- 
ferent non-thermal fluxes. Under the assumptions we made 
the asymmetries that might arise in runaway bowshocks do 
not produce a difference in the emitted spectrum. The ob- 
tained SEDs depend essentially on the particular assump- 
tions made for the particle acceleration, the magnetic field, 
and the dust emission. 

The emission might be detectable at several wave- 
lengths, provided that the source is close enough and long 
exposure times are used - a good candidate is Oph -. The 
synchrotron emission expected at radio wavelengths might 
be detectable, as in the case of BD +43° 3654. The unde- 
tectability can establish constraints on parameters such as 
the magnetic field in the shocked wind. Stellar bowshocks 
might also be detectable at X-ray wavelengths, although no 
runaway bowshock has been observed at these energies so 
far. Finally, a system like C Oph might be detectable at 7- 
rays by the future ground-based detector CTA, as well as by 
the Fermi satellite. The energy range between 1 GeV-1 TeV 
offers the best prospects for the study of runaway stars as 
non-thermal emitters. Our work suggests that bowshocks of 
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Fig. 11. Computed SED for C Oph bowshock, at d - 222 pc. The sensitivity for CTA, Fermi, MAGIC, XMM-Newton 
and VLA. IRAS data are also shown. 



runaway stars might constitute a new class of high-energy 
sources to be explored in the near future. 
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